Norm-conserving Hartree-Fock pseudopotentials and their asymptotic behaviour 
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We investigate tlie properties of norm-conserving pseudopotentials (effective core potentials) gen- 
erated by inversion of the Hartree-Fock equations. In particular we investigate the asymptotic 
behaviour as r — > oo and find that such pseudopotentials are non-local over all space, apart from a 
few special special cases such H and He. Such extreme non-locality leads to a lack of transferability 
and, within periodic boundary conditions, an undefined total energy. The extreme non-locality must 
therefore be removed, and we argue that the best way to accomplish this is a minor relaxation of 
the norm-conservation condition. This is implemented, and pseudopotentials for the atoms H— Ar 
are constructed and tested. 
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The pseudopotential approximation is a vital part of 
the practical application of ah initio methods to problems 
in quantum chemistry and condensed matter physics. In 
the pseudopotential approach the tightly bound core elec- 
trons are removed and their influence on the rest of the 
system is represented by an effective potential. This re- 
duces the number of electrons that must be considered 
and yields a smoother potential. These properties sig- 
nificantly reduce the computational effort required for 
complex systems. 

Over the past few years ultra-soft^'^ pseudopoten- 
tials and separable norm-conserving Kleinman-Bylander"^ 
pseudopotentials have come to dominate within plane- 
wave Density Functional Theory (DFT) technology. 
Within the quantum chemistry community it is more 
usual to use non-separated norm-conserving pseudopo- 
tentials, which are the kind we deal with here. Such 
pseudopotentials can, however, be converted to the sep- 
arable Kleinman-Bylander form. 

Our main interest is in diffusion quantum Monte Carlo 
(DMC) calculations^. Norm-conserving pseudopoten- 
tials are routinely used within DMC calculations^"*', and 
there is evidence to show that Hartree-Fock (HF) theory 
provides better pseudopotentials for this purpose than 
DFT—. HF pseudopotentials have been used in various 
correlated valence quantum chemistry calculations, and 
we suggest they would be suitable for use in perturbation 
theories such as the GW scheme, which has been applied 
to a number of condensed matter systems^. 

It turns out that most of the HF pseudopotentials 
available in the quantum chemistry literature diverge at 
the origin, normally like l/r^ or l/r, which makes them 
unsuitable for use in DMC calculation o^' ' . Re- 
cently Ovcharenko et al}'^ developed HF pseudopoten- 
tials for Be-Ne and Al-Ar which are finite at the origin, 
and we will compare the pseudopotentials generated here 
with these. We would like our pseudopotentials to be as 
smooth as possible, since lack of smoothness can reduce 
the transferability of a pseudopotential. In addition we 
wish to generate HF orbitals using a variety of basis sets, 
including plane waves, for which a smooth pseudopoten- 
tial reduces the computational cost. 



It is advantageous to make the region over which the 
pseudopotential is non-local as small as possible as the 
evaluation of the non-local energy within DMC is expen- 
sive. Having a small region of non-locality also tends to 
promote transferability, but tends to make the pseudopo- 
tential less smooth. 

An important issue arises in generating HF pseudopo- 
tentials by inversion of the HF equations. It has long 
been understood that if the norm is not conserved the 
Hartree potential due to the pseudo-ion does not decay 
as —Z^ai/r at large r (where Zyai is the ionic charge)^. 
It has also been appreciated that the exchange inter- 
action gives a long ranged tail to the pseudopotential 
which should be removed'^^-M. We find that even if norm- 
conservation is enforced the exchange interaction results 
in a "non-Coulombic tail" that decays more slowly than 
the Coulomb part of the pseudopotential, persists far 
from the atom, and is non-local. Non-Coulombic tails in 
HF pseudopotentials have been reported before in the lit- 
erature, but the form that we derive here is different from 
that described by Kahn et alr^ and Hay et ali^, who sug- 
gest that the deviation from —Zyai/r decays faster than 
l/r. 

These non-Coulombic tails increase the cost of calcu- 
lations, lead to unphysical results and, in the case of pe- 
riodic boundary conditions, leave the total energy unde- 
fined. Consequently the tail must be removed for practi- 
cal applications. Although past workers have constructed 
HF pseudopotentials with this long range non-local tail 
removed, we are not aware of a detailed discussion in the 
literature of the form and magnitude of this tail. Un- 
derstanding this long ranged effect is important when 
considering how to remove it. 

In what follows wc investigate the general properties of 
norm-conserving^^ HF pseudopotentials, with particular 
emphasis on their non-locality and asymptotic behaviour. 
We present a method to localise the potential such that 
it is non-local only within a relatively small region sur- 
rounding the nucleus, and enforce the desired asymptotic 
behaviour. 

In section H] we discuss pseudopotentials obtained from 
inversion of the HF equations, in section |TT] we derive 
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the asymptotic behaviour of these pseudopotentials, and 
report detailed results for Ne. In section Hill we present 
and analyse results for the atoms H-Ar. We draw our 
conclusions in section [TVl 

Atomic units are used throughout, unless otherwise in- 
dicated. 



I. PSEUDOPOTENTIALS FROM INVERSION 
OF THE HF EQUATIONS 

The pseudopotential generation procedure we apply is 
similar to that used by others, such as Hay and Wad^, 
or TrouUier and Martins^. In this section we provide 
a brief summary of the procedure for generating norm- 
conserving pseudopotentials, primarily to define the no- 
tation and context. 

For an isolated atom a HF orbitali^, t/jum, may be 
labelled by the quantum numbers i, U and m. Separating 
the angular and radial coordinates leads to 



r 

The radial part of the orbital, (t>i{r)^ satisfies 

h{k + 1) , yeff 



2r2 



where the effective potential for orbital i is 
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z 
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Vh[p] 
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(2) 



(3) 



The first term arises from the nuclear charge and the 
second is the Hartree potential due to the total electron 
density. The third term is the effective exchange poten- 
tial, which is different for each orbital. For convenience 
in what follows we have included the (cancelling) self- 
interaction in the second and third terms. 

To construct a pseudopotential we first partition the 
atomic orbitals into core orbitals, whose influence will 
be represented by the pseudopotential, and valence or- 
bitals, which will be represented by the pseudo-orbitals. 
The eigenvalues of the pseudo-orbitals are constrained to 
equal the equivalent all-electron (AE) eigenvalues, and 
the pseudo-orbitals themselves are constrained to equal 
the equivalent AE orbitals outside of a "core radius", rd- 
Inside rd the orbitals are given by an as-yet unspecified 
function /i, so that 



Mr) 



_ / Mr) r <rc 
(j)i[r) r >rc 



(4) 



fi{r) is chosen such that (pi{r) is nodeless, smooth up to 
a certain order of differentiation at r = r^i, and that its 
norm is equal to that of 0i(r). 

The pseudo-orbital for state i satisfies 



1 kik + 1) yeff 

2 dr2 2r2 ' 



(5) 



and hence the effective potential, V^^^ , can be obtained 
by inverting Eq. ([5]). The pseudopotential for state i, 
Vi(r), is then defined by the pseudo-atom equivalent of 
Eq. ©: 



V,{r)+Vn[p\ 



(6) 



To use the pseudopotential in a calculation for a 
molecule or solid it is expressed in terms of projection 
operators, and separated into local and non-local parts. 



^pseudo — ^local{r) 

+ flY. \yim)iVl{r)-Vlocalir)){Yir, 



I 7n— — l 



(7) 



The orbitals with I > Imax feel the local potential, Viocai ■ 
Here the index I is interchangeable with the pseudo-atom 
orbital index, i, since the valence states chosen to con- 
struct the pseudopotential must have unique / quantum 
numbers. 



II. EXTREME NON-LOCALITY OF HF 
PSEUDOPOTENTIALS 

In this section we address the question of the locality 
of the pseudopotentials defined above. We define a pseu- 
dopotential to be local within a radius rio^ if, for all i 
and j, 



\V^{r)~VJ{r)\ <6iorr>rio 



(8) 



for some rioc and small S. 

By definition 0j ~ cj)i for r > and hence V.^^'^'^ — 
This equality, together with Eqs. ([3]) 



,eff 
i 

and ([6]) gives 



Vr{r) 



Z 



Vh[p-~p] 



0, 



r > Tn 



(9) 



The first two terms are independent of the orbital and, 
since p and p differ by the exponentially decaying core 
contribution, the sum of these terms approaches —Zyai/r 
as r — + oo. However, since the HF exchange potential is 
non-local there is no reason to expect the final two terms 
to cancel, or to approach zero faster than 1/r as r — s- oo. 

To obtain the asymptotic form of t4[{(/'}, ^i](/'i/</'i we 
start from the asymptotic behaviour of the atomic HF 
orbitals. Handy et al^ have shown that the asymptotic 
form of the atomic HF orbitals as r oo is dominated 
by an exponential decay with exponent a = {—ieuo)^ i 
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where HO denotes the highest occupied orbital. Handler 
et al?^ derived the full asymptotic form, 

(t>r{r) - rP'+\a, + b,r-^ + . . .)e-"^ (10) 

where we use ~ to denote "asymptotically approaches as 
r oo". They found expressions for f3i which can be 
summarised as 

J k §1 

= HO =Iho /3 

y^HO ^Iho P~\k-lHo\-l (11) 
^HO =Iho^O /3 - 3 

y^HO =IhO=0 13- 2{lrmn + 1) 

where f3 — {Z — {2N — 1) — a) / a, there are 27V electrons in 
the atomic configuration, and Imin is the lowest non-zero 
I value in the configuration. 

In Eq. pT|) and the rest of this section we limit our 
analysis to closed shell atoms to aid clarity. The gen- 
eralisation to open shell atoms using spin averaging or 
eigenfunctions of the total angular momentum is straight- 
forward but algebraically complex, and does not change 
our conclusions. 

Eqs. (fTO)l and (fTTj) do not apply to s-electron-only 
atoms. In this case the exponent for each orbital is dif- 



ferent and is given by = (— 2ei)^, and it is straightfor- 
ward to demonstrate that the sum of the exchange terms 
in Eq. ([9|) approaches zero exponentially as r increases. 
We do not consider these pure-s atoms in what follows. 

Expressing the HE exchange potential in spherical po- 
lar co-ordinates and taking the limit as r — > oo gives 

V,m,k]Mr) - -^5](2/, + l)u;, (^J; ^ 

where the Wigner 3j notation has been used for triple 
integrals over spherical harmonics. The occupation of 
shell j is denoted by Wj , where Wj = 1 for a full shell. 
We define the effective exchange potential for orbital i 

as 

W = —r\ ■ 13) 

Next we explicitly separate the exchange interaction be- 
tween electrons within shell i from the remaining ex- 
change, and use the asymptotic forms for the orbitals 
of Eq. Uni), to obtain 



where riij = Pj — f3i — I — 1. The first term in Eq. is 
always present and is a part of the self-interaction correc- 
tion (SIC) which is the same for all orbitals. The second 
term is part SIC and part exchange interaction, and is 
non-zero only for li > 0. If present this introduces terms 
cx i/r^/c-i-i ^Qj. ^ _ 1, . . . , which are different for each 
orbital. 

where m — minj^^/^^ ~ Iho\ and Imin is as defined for 

Eq. dm). 

Using Eq. we may write the effective exchange 
potential in the limit r — s- cx) as 

The third term is the remainder of the exchange inter- 
action, and we consider the maximum power riij that 
appears in this expression with a non-zero coefficient. 
Given that the Wigner 3j symbol is non-zero only for 
\li — lj\<l<li+ Ij and li + lj+l even, and that the inte- 
gral is zero for 1 = 0, we find that the maximum riij = i^i 



takes the values 

= H0 =Iho -(2m + 2) m 
^ HO ^Iho \k-lHo\ 

(15) 

^HO =Iho^O 2 

7^ HO = IhO =0 Imin 



r'-4>j(j3idr 



{ij}\n,,=o 

poo 



(16) 
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where the Wigner 3j and its pre-factor in Eq. ((T4)l have 
been subsumed into the parameters a^i and biji, and the 
sums are taken over the I and j that give Uij = and 



= — 1. For i = HO, only the SIC term appears since 



= h 



'ijl 



0. Note that for i ^ HO, V^{oo) is non-zero 



and the effective exchange potential contains a term 
in addition to the SIC. 

Using the same line of argument, the pseudo-atom ef- 
fective exchange potential at large distances is given by 



1 /■°° 
{ij}|„,,=o •'° 



(17) 



again with Oy; — biji = for i — HO. The sum over l,j 
is as before, and a similar non-zero limit and term 
appear. The occupation of pseudo-shell j, Wj, is defined 
as for the AE atom, but is zero for core states. Since 
the integrals in Eq. and Eq. pT)) are not in general 
equal, the limits for the AE and pseudo-atom effective 
exchange potentials are not equal and do not cancel in 
Eq. In addition, exchange terms due to core- valence 
interactions will not cancel as no core orbitals are present 
in the pseudo-atom. 

The asymptotic behaviour of the pseudopotential can 
be obtained to 0(r^^) from the explicit forms of the 
external, Hartree, and exchange potentials. Starting 
with the AE effective potential, Eq. ^ can be expressed 



{Z - {2N- l))r- 



r (j)j(j)idr 



|n„=-l 



{lj}\nij = 



(18) 



If we write the asymptotic form of the pseudopotential 
as 



V, ~ s, - [Z^ai + g^)r-' + 0(r-'), (19) 

where s, and gi are constants then the equivalent expres- 
sion for the pseudo-atom effective potential is 

V!^^ ~ s,-(Z™j-(27V,,,-l)+.g,)r-i 



{l]}\n,j=0 



{lj}\nij=-l 



Wjaijir ^ / r'-(f>j(f>idr + 0{r ^) , 



(20) 



where there are 2Nyai valence electrons. Since 
y^ff for r > we may equate the coefhcients of 
and r° to give 



= HO 



E{y}|„.,=oay' Wj r'4>j4),dr - Wj J^^ r^(l)j(j),dr i ^ HO , 



(21) 



and, assuming that the charge of the AE and pseudo-atoms are equal {Z — 2N ~ Zy^i — "^N^ai), 

i:^ HO 

9i = 



Wj /o°° r^(j)j(t)idr - wj r'-(j)j(j),dr 
I 



(22) 



i 7^ HO 



For i 7^ HO, Si and gi are generally non-zero since the 
AE and pseudo-orbitals differ in the core region. In ad- 
dition, terms due to exchange between state i and core 
states have no counterpart in the pseudo-atom to cancel 
with, as Wj — 1 and wj = for the core states. This 
orbital dependent "offset potential" and "ghost charge" 
result from a remnant of the AE exchange interaction 
remaining in the pseudopotential as it is not completely 
removed by the exchange interaction in the pseudo-atom, 
and is generally due to both valence- valence and valence- 



core exchange. Because these terms are orbital depen- 
dent the pseudopotential itself is non-local over all space 
and does not satisfy the locality criteria of Eq. ^ for any 
rioc- It should be stressed that this extreme non- locality 
is necessary if we require our HF pseudopotential to be 
both norm conserving and to reproduce the AE valence 
eigenvalues and orbitals outside of the core region. 

To illustrate the physical consequences of this be- 
haviour we consider the total energy of a collection of 
identical pseudo-atoms. The total energy should tend 
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towards the total energy of the isolated pseudo-atoms as 
they are separated further apart. This is true if each 
channel of the pseudopotential approaches —Z-uai/r for 
large r. For the extreme non-local behaviour found here 
this is not the case, and an unphysical interaction persists 
at large distances. For extended systems this is catas- 
trophic because the total energy is undefined. 

Extreme non-locality may be avoided in the construc- 
tion of a pseudopotential by using a different configura- 
tion to construct each channel, such that the configura- 
tion used to construct channel I has I = I ho- Although 
this is practical for H and He, for heavier atoms this 
requires the construction of pseudopotentials from elec- 
tronic configurations far from those we wish to use them 
for. It would not be reasonable to expect the resulting 
pseudopotentials to be transferable to systems of inter- 
est, and they would not reproduce the valence states of 
the neutral isolated atom. 



A. The Kleinman-By lander form 

The Kleinman-Bylander form is expressed in terms of 
operators that project the orbitals onto a set of basis 
functions, and takes the general form 

V'''' = Viocai{r) + |^,;)A,(^,| , (23) 

for some {'(/'}• The pseudopotentials discussed in section 
H] may be expressed in this form if we define 

D^J = {^^\VJ-Vlocal\^J), (24) 

where the ipi are the eigenstates of the pseudo-atom. If 
a complete set of V'i is included in the sum in Eq. (|23p . 
this form would be exactly equal to Eq. fT]). In prac- 
tice a more useful approach is to use a small number 
of localised, bound, atomic states (such as those used 
to construct the original pseudopotential) with Viocai{f) 
chosen such that Viocai —Zyai/r as r —>■ oo. For this 
choice of localised bound states, V^^ is well-behaved at 
large r, since the subspace covered by the projection does 
not include this asymptotic region - the non-local part of 
the pseudopotential has been localised by representing it 
within a locahsed subspace. 

We conclude that converting a norm-conserving HF 
pseudopotential to Kleinman-Bylander form removes the 
long ranged non-locality, and so allows these potentials 
to be implemented computationally in a manner which 
avoids the unphysical behaviour described previously. It 
should be noted that the "frozen remnant" of the atomic 
exchange interaction which makes up the non-Coulombic 
tail is still present (although limited to a subspace rep- 
resented by the pseudo-orbitals), and it may still lead to 
significant errors. 



B. DFT pseudopotentials 

Extreme non-locality does not occur for pseudopoten- 
tials constructed within Kohn-Sham DFT because the 
effective potential that represents the electron-electron 
interaction is the same for all electrons (i.e., it is local). 
More explicitly, in DFT the local exchange-correlation 
potential takes the place of the exchange potential in 
Eq. ^ and hence the equivalent expression for the pseu- 
dopotential is 

Vdr) = -- + Vh[p-p] + V^c[p]-V^e[p], r>r„. (25) 
r 

For r > Tc, where Vc — max [rd], all Vi{r) are the same, 
so the pseudopotential is local outside of the core region. 

Within the Optimised Potential Method {OFM)^^ 
the exchange-correlation energy functional is replaced by 
an "exact" exchange energy functional (possibly with an 
added approximate correlation functional) defined explic- 
itly in terms of the Kohn-Sham orbitals and implicitly in 
terms of the charge density. This exact exchange (EXX) 
functional takes the form of the Fock term of HF, but is 
evaluated using the KS orbitals and hence the calculation 
remains KS-DFT^** and the 0PM exchange energy does 
not equal the HF exchange energy. Norm-conserving 
pseudopotentials have been defined within this theoret- 
ical picture by Engel et al. and others, and "spurious 
long range structure"— i^i^l consistently appears in these 
pseudopotentials . 

Given the highly non-local nature of the Fock term it 
seems reasonable that the exchange potential is also an 
extremely non-local functional of the density. We also ex- 
pect the exchange potential to approach zero as r — > oo, 
since this behaviour is part of the definition of the EXX 
potential itself^^. These two observations, together with 
Eq. ([25| . explain the existence of the structure found 
by Engel et al.'^'L They found that pseudopotentials con- 
structed by removing this long range structure resulted in 
improved bond distances and energies for simple dimers 
and bulk Al. Note that the 0PM case is in some regards 
similar to that found in the HF case, but it is essentially 
different in that the potential is local for r > Vc and 
approaches —Zyai/r faster than 1/r. 

For local (e.g., the local density approximation or 
LDA) or semi-local (e.g., the generalised gradient ap- 
proximation or GGA) exchange-correlation functionals 
we note that p and p differ by the exponentially decaying 
core contribution and obtain the relation 

lim V,{r) = -— + 0{e-"^). (26) 

The absence of an offset or ghost charge is a consequence 
of the locality of the approximate exchange-correlation 
functional. It is clear that there is no DFT analogue of 
the extreme non-locality present in HF pseudopotentials. 
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III. RESULTS 

A. Size of the extreme non-local terms 

To construct a pseudopotential from an AE atomic 
HF calculation we must define a form for the pseudo- 
orbitals inside the core radius, fi. We chose the TrouUier- 
Martinsi^ form 



Mr) 



exp 



m— 



(27) 



which is based on the observation that the pseudo- 
orbitals are smoother if the odd derivatives of the 
screened pseudopotential at the nucleus are set to zero. 
The coefficients are determined by requiring norm- 
conservation of the charge in the core region, continuity 
of and its first four derivatives at Vd, and the require- 
ment that Vi'^^^ has zero curvature at the origin. 

As an example we consider Ne (ls^2s^2p^) and perform 
a numerical HF calculation and construct (j) for the 2s and 
2p orbitals. We then obtain s and p pseudopotentials as 
described in section |T1 We also calculate the r~-^ and r° 
contributions to the asymptotic behaviour directly from 
the pseudo-orbitals. 

First we note that, for Ne, the highest occupied orbital 
is i = 2p. From Eqs. ([14|) and Eq. HSl) it follows that 
the p part of the pseudopotential approaches ~Zyai/r as 
due to the 2p — 2p exchange interaction. 

For the s part, the exchange interaction contributes 
the 2s — 2p term, which is the only interaction appearing 
in Eqs. (|18p and ([20]). To obtain values for S2s and g2s we 
first obtained the asymptotic behaviour of the effective 
potential V^'^'^{r) by direct substitution of the asymp- 
totic form of the orbital into Eq. ([2|). Equating powers 
of r leads to 



(e2s - e2p) = -a2p,2s,i / r4)2p<p2sdr , (28) 
Jo 



and 



2V2\e2p\^ ^ b2p,2s,i r(t)2p(t)2sdr . (29) 

Using the numerical values of the integrals and eigen- 
values we obtain the values of a2p,2s,i and &2p,2s,i, and 
by using Eqs. dH]) and (QH) we obtain S2s = 0.0095 
a.u. (0.26 eV) and g2s — 0.023 e, that is, a positive offset 
to the potential and a positive ghost charge. 

The difference between the s part of the Ne pseudopo- 
tential and the Coulomb potential, Vs -\- Z^ai /r is shown 
in Fig. [TJ It is clear that the pseudopotential does not 
approach the Coulomb potential for large r. The calcu- 
lated asymptotic form, S2s-\-g2s/r, is also shown. The dif- 
ference between the pseudopotential and the calculated 
asymptotic form at large r arises mainly from a 1 /r^ term 
which was not included in the asymptotic form. 
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FIG. 1: The difference between the s part of the Ne pseudopo- 
tential and the Coulomb potential, Vs + Zvai/r, against ra- 
dial distance (solid line), and the calculated asymptotic form 
(dashed line). 



The magnitude of the deviation from the Coulomb po- 
tential at large r shown in Fig. [1] is typical of that found 
for all p and d-block atoms. We conclude that the non- 
locality persists far from the nucleus and is of signifi- 
cant magnitude for these atoms. For s-block atoms the s 
channel is due to the HO orbital, hence sho = Qho — 0- 
The p and d channels are generally obtained from excited 
atomic states for which the offset is typically ^ I0~^ a.u. 
for the channel which does not correspond to the HO 
orbital. 



B. Localising the pseudopotential 

Here we consider several methods for removing the 
non-Coulombic tails of the pseudopotentials in order to 
localise them within some radius rioc- In what follows 
we define the localisation radius, rjoc, as in Eq. ([H]), with 
,5 = 5 X 10^^ a.u. 

Previous workers (for example. Greet and Lester-i) have 
expanded pseudopotentials in a Gaussian basis, a stan- 
dard form for quantum chemistry applications. Due to 
the functional form of this expansion the parameterised 
pseudopotential is non-local in a finite region surround- 
ing the nucleus, and approaches the Coulombic form as 
r ^ oo, hence the expansion may be used as a locali- 
sation procedure. Generally this expansion is obtained 
numerically by requiring the parameterised pseudopo- 
tential to reproduce the pseudo-states of the original 
pseudopotential to within a given accuracyi^ Although 
this method can work well it is difficult to control the 
smoothness and the localisation radius, r/oc, since the 
error in the expanded pseudopotential (and accompany- 
ing pseudo-states) must compensate for the absence of 



7 



the long range behaviour present in the original pseu- 
dopotential. We prefer to control the localisation of the 
original pseudopotential and then, if desired, expand the 
resulting potential in a Gaussian basis. Gaussian expan- 
sions are not considered further in this paper, but we will 
present parameterised pseudopotentials obtained in this 
way in a separate paper^^. 

We have freedom in our choice of Td and fi (r) , and by 
varying these quantities we might be able to localise the 
pseudopotential in the sense of Eq. ([5]) , since variations 
in these quantities will alter the integrals in Eqs. ((2T|) 
and ([22]) . For Ne (and for most atoms) the first term in 
Eqs. (PT|) and is larger than the second because the 
AE orbitals possess radial intervals where they are neg- 
ative, while the pseudo-orbitals are positive for all r. In 
order to decrease the value of this integral, causing S2s 
and g2s to become closer to zero, fi{r) might be chosen 
such that, on average, the electron is closer to the nu- 
cleus. This might be achieved by increasing rd, allowing 
the maximum of 4>2s to be closer to the nucleus. Also, the 
form of fi (r) may be altered (a convenient way of achiev- 
ing this is to relax the continuity of the 3'''' and/or 4*^* 
derivatives). Calculations suggest that neither of these 
approaches is useful - for a good cancellation of the inte- 
grals in Eqs. (|2H22p rd must be very large, or fi{r) must 
take a form such that the corresponding pseudopotential 
is far from smooth. 

Next we consider an elegant procedure described by 
Engel et alM and used by them to remove the (local) 
intermediate range structure present in 0PM pseudopo- 
tentials. This approach employs a self-consistent modi- 
fication of the pseudo-orbitals over all space, while con- 
serving the norm. The same method may be applied to 
localise HP pseudopotentials. 

The pseudopotential at iteration n is denoted by Vi-m 
and the associated pseudo-orbitals are {4>i}n- The pseu- 
dopotentials calculated in section |T] correspond to n = 1. 
We construct an effective potential which is not yet self- 
consistent, 

y^■n = yVh[Pcor,\+Vh\PnU 5 , (30) 

where —Z/r + Vh[pcore] is the Coulomb+Hartree poten- 
tial due to the ionic core. Equation ([50]) is the external 
potential we desire outside of the core radius together 
with the effective interaction potential of the pseudo- 
orbitals {(f>i}n- Taking the correct asymptotic form for 
the orbitals, we integrate in from oo to rd at energy et to 
obtain the new orbitals {4)i}n+i in this region, and these 
are normalised by requiring that the norm of the newly 
generated and original pseudo-orbitals outside of rd are 
equal. In the core region the new orbital is constructed by 
using the TrouUier-Martins form and criteria described at 
the beginning of this section. Taking {0i}„+i, we invert 
the HP equations as described in section[T]to find the new 
external potential that results in these orbitals {0i}„+i 
on solution of the HP equations, Vi-_n+i- 



The new set of orbitals are then used in Eq. ((30)) to 
generate the effective potential V^[^i and the entire pro- 
cess repeated until self-consistency is achieved - we use 
|Vi;ri+i — Vi-n\ < 10~^ a.u. for all r. At self-consistency 
Vi-n is the pseudopotential we require. This method finds 
the HP orbitals that take the TrouUier-Martins form in 
the core region, result from a local ionic potential out- 
side of the core region, and have the same eigenvalues 
and norm as the original AE atom valence states. 

A number of difficulties presented themselves with this 
method. We found that for atoms with few valence elec- 
trons it was not possible to achieve self-consistency. This 
appeared to be because the analytic form used for the or- 
bitals in the core region was not general enough to give 
a self-consistent solution. Furthermore, although Engel's 
approach conserves the norm of the original AE states, it 
does not conserve the magnitude or the radial derivative 
of the orbitals at rd- Consequently the new pseudopo- 
tential does not conserve the logarithmic derivative or 
its energy derivative, and does not reproduce the scat- 
tering properties of the AE atom. For example, for the 
s-channel of Ne the new pseudopotentials results in a 
—0.24% change in the magnitude of the orbital at the 
core radius, a -f3.9% change in the logarithmic deriva- 
tive, and a -1-0.40% change in the energy derivative of the 
logarithmic derivative. These errors may be acceptable, 
although the error in the logarithmic derivative might be 
problematic. 

By introducing an interface region between the core 
region and the valence region (where the pseudopoten- 
tial is ionic) it is in principle possible to modify Engel's 
method such that the magnitude of the orbital, the log- 
arithmic derivative and the energy derivative of the log- 
arithmic derivative are all conserved at the core radius. 
This was attempted using the TrouUier-Martins form to- 
gether with additional criteria. We found that either self- 
consistency could not be achieved or that the resulting 
pseudopotentials were far from smooth. 

The method we chose to use is to apply a transforma- 
tion to the original pseudopotential such that deviations 
from the Hartree-|-Coulombic ionic core potential are re- 
moved beyond a certain radius but desirable properties 
of the original potentials are preserved. We investigated 
the transformation 

'%{r) + V,{r) r<R, 

X [it (r) + v., (r) - Vh [pcore] + f ) 
. + Vh [pcore] - f r> Ri , 

(31) 

where — Z / r + Vh[pcore] is the ionic potential. The Gaus- 
sian decay of the non-Coulombic part of the potential 
occurs over a distance ~ l/a~2 , and the zeroth and first 
derivatives of V^'°'^ are continuous at Ri. The precise 
value of a has very little influence on the final result - 
the cutoff radius, governs the quality of the pseu- 
dopotential. 
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For the simple cutoff function corresponding to 7; = 
and a~2 = 0.15 a.u. we chose Ri as the minimum value 
for which the je^ — < 10^^ a.u., where the e^°^ are 
the eigenvalues of the new pseudopotential. For Ne this 
new pseudopotential is local (in the sense of Eq. ([8])) 
outside of a radius of rioc — 3.35 a.u. Although this is 
large in comparison with the core radii of the orbitals 
themselves, it may be acceptable, although it would be 
expensive for DMC. Moreover, this approach enforces a 
Gaussian decay on a monotonically increasing function, 
essentially creating a (smoothed) step-like behaviour in 
the pseudopotential, and it seems likely that this would 
introduce transferability problems. Note that if we seek 
to remove this step-like behaviour by decreasing a then 
rioc rapidly becomes impractically large. 

A more successful approach is to choose Ri = rd and 
— rc/16 (as before, the final result is insensitive to 
the choice of this parameter) so that the pseudopotential 
is non-local in a small region and no step-like behaviour 
occurs. In order to reproduce the eigenvalues of the orig- 
inal pseudopotential, ji must be non-zero, and we use 



r<R^ 
r>R^, 



(32) 



where a is chosen as before and qi and pi are orbital 
dependent parameters. This function is continuous and 
smooth at r — Ri, and the second derivative is zero at 
r = (see the discussion following Eq. ([?f|) ). Our goal 
is to search for values of these parameters such that the 
localised pseudopotential preserves certain desirable fea- 
tures of the original pseudopotential. We require that 
the conditions 



or 



Joe 



dr 



~lo 



= 0, 

- = , 



are satisfied in order to reproduce the original eigenvalues 
and aid transferability. Conservation of these two quan- 
tities takes precedence over norm conservation and we do 
not require the norm of ^'"'^ to be equal to the norm of 
(j)i as we could not find a reliable method to achieve this 
and give a smooth pseudopotential. 

Conservation of the logarithmic derivative of each 
pseudo-orbital at Ri may be achieved exactly, but a sim- 
pler and more convenient approximation is used here. 
Starting with the equations for the original pseudo- 
orbital and the pseudo-orbital resulting from the localised 
pseudopotential we may follow a similar derivation to 
that used to arrive at the standard norm-conservation 
condition to obtain 



— In (^'"'^ - — In ( 
dr ' dr 



R^ 



VLoc 



AV(l)idr - {e: 



-loc 



lloc [, 



idr 



(34) 
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FIG. 2: The localised pseudopotential for Ne. The inset shows 
the difference between the localised and original pseudopoten- 
tials. 



where the eigenvalues are not yet constrained to be equal 
and Ay = y^^f^^"'^ _ V^^^ is the difference between the 
effective potentials for calculations carried out with the 
original and localised pseudopotentials (for the configu- 
ration used to construct the pseudopotentials). Our ap- 
proximation is to take AF to be close to the difference 
between the external potentials. A]/ « 1//°^ — Vi = 7i, 
for r < Ri. 

We start by using simple iteration. For a fixed qi 
we take an initial value of = and solve for the 
resulting pseudo-orbitals. We then update pi by solv- 
ing J^' 4>'f'^ji(f>^dr — for pi and continue until self- 
consistency is achieved. This is repeated for different 



(33) values of qi and the condition e: 



~loc 



= is enforced 
by straightforward bisection. This process is repeated 
to obtain self-consistency across all channels, and we 
obtain eigenvalues which differ by < 10^^ a.u., and 
/o^' 4'^i"^li4'idr = to numerical precision. Note that as 
we have not taken into account the change in the effective 
potential due to changes in the exchange and Hartree po- 
tentials, the LHS of Eq. ([M)) is not zero, and there will 
be a small change in the logarithmic derivative. 

Localised pseudopotentials for Ne are shown in Fig. [51 
The s and p channels were obtained from the neutral 
ground state, while the d part was obtained from an ex- 
cited ionic state (Is22s^2p2-^^3d°-^^). All excited state 
configurations used to obtain pseudopotential channels 
not bound in the ground state configuration are taken 
from Bachelet et al^ Inset in the same figure is the dif- 
ference between the localised pseudopotential and the 
original numerically exact pseudopotential. This figure 
demonstrates that our localisation procedure deforms the 
pseudopotentials by a small amount, mostly within the 
core region. 

For all atoms considered we found that this process 
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FIG. 3: Difference between the charge density resulting from 
the localised (p'°'^) and original (p) pseudopotentials for Ne. 
To achieve localisation a small amount of charge near the 
nucleus is redistributed. 



results in a localised pseudopotential which modifies the 
form of the associated orbitals over all space, but by a 
very small amount - about 0.1% of each electron is re- 
moved from a region around the core radius and redis- 
tributed into two regions, one closer to and one further 
from the nucleus. The removal of offset and ghost charge 
(and higher order terms) results from a small change in 
the asymptotic behaviour of the orbital at large r. Nor- 
malisation must be preserved, and hence the orbitals are 
also changed slightly at small r. This effect is illustrated 
in Fig. [31 which shows the difference between the radial 
charge density of Ne resulting from the localised pseu- 
dopotential and the radial charge density resulting from 
the original pseudopotential. 

The transformed Ne pseudopotential is local (in the 
sense of Eq. ([8])) outside of a radius of rjoc = 0.74 a.u. 
A similarly localised pseudopotential constructed by ap- 
plying Eq. ([3T|l with ji — (using the same a and Ri 
values) results in a large error in the 2s eigenvalue of 
-1.2 X 10~2 a.u. 

As the localised version of the pseudopotential is not 
equal to the original pseudopotential it is not exactly 
norm conserving, and does not preserve the form of or- 
bitals outside the core region or at the core radius. To 
illustrate the magnitude of these effects we return to Ne, 
specifically the s channel as this is the major source of the 
extreme non-locality present in the original pseudopo- 
tential. The localisation process results in a —0.27% 
change in the magnitude of the orbital at the core ra- 
dius, a —0.05% change in the logarithmic derivative (due 
to assuming that the exchange and Hartree terms to do 
not contribute in Eq. ([M)) ). and a -1-0.80% change in the 
energy derivative of the logarithmic derivative. 

The minor relaxation of norm-conservation that ac- 



Atom 






Ted 


rioc 


H 


0.50 


0.50 


0.50 


0.49 


He 


0.60 


0.60 


0.60 


0.59 


Li 


2.19 


2.37 


2.37 


2.71 


Be 


1.88 


1.96 


1.96 


2.20 


B 


1.41 


1.41 


1.41 


1.63 


C 


1.10 


1.10 


1.10 


1.29 


N 


0.94 


0.88 


0.84 


1.09 





0.80 


0.75 


0.99 


1.13 


F 


0.70 


0.64 


0.89 


1.02 


Ne 


0.63 


0.57 


0.63 


0.74 


Na 


2.70 


2.85 


2.85 


3.25 


Mg 


2.38 


2.38 


2.38 


2.71 


Al 


1.94 


2.28 


2.28 


2.64 


Si 


1.67 


2.01 


2.06 


2.36 


P 


1.48 


1.71 


1.71 


1.98 


S 


1.33 


1.50 


1.50 


1.74 


01 


1.19 


1.34 


1.34 


1.55 


Ar 


1.09 


1.20 


1.31 


1.54 



TABLE L Core radii and localisation radii (a.u. 



companies the localisation may influence the transferabil- 
ity of the pseudopotential. However, the familiar conser- 
vation of the first order change in the radial logarithmic 
derivative with respect to energy is only valid when the 
potential inside the core region is fixed. ^° For both the 
HF and DFT this is not the case even for an isolated 
atom, since V^^^ is a functional of {0} and each (pi is a 
function of {e}, and we have shown (from Eq. (|34p ) that 
the norm-conservation relation is replaced by 



d_d_ 

dej dr 



h\\(j)i 



2 



\<i>,\^dr 



\<i)i{rciW Jo 



' dr. 



(35) 



This applies to both the AE and pseudo-atoms. Given 
that Eq. (|35|l tells us that first order changes in the ra- 
dial logarithmic derivative depend on the variation of the 
effective potential with as well as the norm, and that 
norm-conserving pseudopotentials have been successful 
in the past, it is likely that the very small relaxation of 
norm-conservation is not significant. 

As a further test of the approach described above we 
have generated localised pseudopotentials for H-Ar. The 
core radii and localisation radii of these pseudopotentials 
are given in TablelH Useful pseudopotentials are required 
to be transferable, that is, to accurately reproduce the be- 
haviour of the AE atom in other environments. In Table 
|TT]we compare excitation energies of the atoms H-Ar re- 
sulting from AE, original pseudopotentials, and localised 
pseudopotential calculations. Excitation energies are cal- 
culated for each atom as the difference between the total 
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energies of the excited and ground states for each type of 
calculation. 

We compare our pseudopotentials with the HF pseu- 
dopotentials of Ovcharenko et al}^. This is the most 
appropriate comparison as both sets of pseudopotentials 
are soft (in the sense of possessing no singularity at the 
origin). The pseudopotentials of Ovcharenko et alM^ are, 
however, parameterised, and we note that some of the 
differences in accuracy will correspond to errors in the 
parameterisation rather than differences in the methods 
of generation. 

In Table IIIII we compare the ionisation potentials 
resulting from the AE, the original pseudopotentials, 
the localised pseudopotentials, and the pseudopotentials 
published by Ovcharenko et al}- in the same manner as 
Table |TT1 Ionisation potentials are defined for each atom 
as the difference between the total energies of the ionised 
and ground states for each type of calculation. 

From Table [Til smd IIIII it is apparent that the differ- 
ences between the energies resulting from the localised 
and original pseudopotentials is negligible in comparison 
to the differences between the energies obtained from the 
AE and pseudopotential calculations. This suggests that 
the localisation procedure has been successful, and that 
the transferability is unlikely to be improved by requiring 
a more exact reproduction of the original pseudo-orbitals. 

Our localised pseudopotentials reproduce the AE ex- 
citation energies well, suggesting that our pseudopoten- 
tials are transferable, and that the small relaxation of 
the norm-conservation criterion has no significant effect. 
Our pseudopotentials also appear to reproduce the exci- 
tation energies rather better than those of Ovcharenko 
et al. In a few cases this appears to be due to the lack 
of a d channel in the pseudopotentials of Ovcharenko et 
al. (e.g. P), but for most atoms considered this is not 
significant (e.g. Al). 

In addition to this our pseudopotentials are softer than 
those of Ovcharenko et ai, see Fig. 2] for a comparison of 
the Si pseudopotentials. 



IV. CONCLUSIONS 

We have shown that in general a non-local norm- 
conserving pseudopotential constructed within HE the- 
ory will be non-local over all space, resulting in a long 
range interaction between atoms that is incorrect and 
not present in the equivalent AE calculation. In particu- 
lar it is apparent that the total-energy cannot be defined 
for extended systems. Although our analysis of extreme 
non-locality was for closed shell atoms within HE theory, 
it is apparent that the same non-locality persists for open 
shell atoms and for Dirac-Fock theory, since the nature 
of the exchange interaction remains the same (although 
the explicit equations become more complex) - it is rep- 
resented by a non-local functional of all states which is 
different for each state. 

This extreme non-locality is of particular relevance 




r (a.u.) 



FIG. 4: Comparison of the s and p channels of the localised 
Si pseudopotentials generated within this paper (Our) and by 
Ovcharenko et al}'^ (OAL). The OAL pseudopotential does 
not have a d channel for us to compare with. 



to the application of pseudopotentials within the DMC 
method. HE pseudopotentials have been shown to give 
good results within DMC calculations, but computa- 
tional efficiency requires the region of non-locality to be 
as small as possible. In light of this we have imple- 
mented a localisation procedure which generates a new 
pseudopotential that is non-local in a small region sur- 
rounding each nucleus. Localisation together with repro- 
duction of the atomic eigenvalues and logarithmic deriva- 
tives is achieved at the cost of a small relaxation of the 
norm-conservation condition, due to an alteration of the 
asymptotic form of the pseudo-orbital. 

We have calculated localised pseudopotentials for H- 
Ar, and tests for these atoms demonstrate the reliabil- 
ity of our approach. We have also demonstrated their 
transferability by presenting calculations of excited states 
for these atoms. This suggests that the approach given 
here provides accurate, well localised, smooth, HF pseu- 
dopotentials suitable for use in ab initio methods such as 
DMC, GW, and quantum chemistry techniques. 

In a future paper we will further generalise our ap- 
proach to relativistic HE pseudopotentials (relativistic 
effective potentials), parameterise these pseudopotential 
in terms of a Gaussian basis set for use in quantum chem- 
istry codes, and present a library of pseudopotentials for 
Z = 1 - 56 and 71 - 80. 
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0.0000 

0.200 ( 
0.4269 


0.0000 

0.4257 


CI 


Ss^Sp^pP] 

Qo2q„4oj1 r4 pi 

Oo Oti [ JT^ J 

3s'3p«pS] 


0.0000 

U.O ( oo 

0.5653 


0.0000 

U.O 1 OtJ 

0.5671 


0.0000 

U.O 1 04: 

0.5671 


0.0000 

U.O ( 04; 

0.5675 


Ar 


Ss^SpSpS] 
3s^3p^3d' pF] 
3s'3p«3d'pD] 


0.0000 
0.4824 

1.1597 


0.0000 
0.4838 

1.1623 


n nnnn 
0.0000 

0.4837 

1.1621 


n nnnn 
0.0000 

0.4855 
1.1643 


Average error 
Maximum error 




+1 X 10"" 
-5 X 10^3 


+7 X 10"^ 
-5 X 10"3 


+9 X 10"'' 
-1 X 10"^ 



TABLE II: Comparison of excitation energies obtained from AE calculations {E'^^), from calculations using the extremely 
non-local and localised pseudopotentials generated within this paper {Eexc and Eix'c respectively) , and from calculations using 
the pseudopotentials generated by Ovcharenko et ali^ {E^^'") (a.u.). 



Atom 




Eton 




^OAL 


H 


0.50000 


0.50000 


0.50000 




He 


0.86168 


0.86230 


0.86230 




Li 


0.19631 


0.19632 


0.19632 




Be 


0.29563 


0.29578 


0.29578 




B 


0.29149 


0.29199 


0.29196 


0.29173 


C 


0.39640 


0.39702 


0.39693 


0.39696 


N 


0.51293 


0.51357 


0.51352 


0.51378 


O 


0.43679 


0.43672 


0.43680 


0.43679 


F 


0.57763 


0.57790 


0.57794 


0.57672 


Ne 


0.72928 


0.72982 


0.72982 


0.73124 


Na 


0.18195 


0.18210 


0.18210 




Mg 


0.24283 


0.24299 


0.24303 




Al 


0.20204 


0.20227 


0.20227 


0.20531 


Si 


0.28123 


0.28175 


0.28172 


0.28164 


P 


0.36901 


0.36984 


0.36977 


0.36957 


S 


0.33171 


0.33061 


0.33060 


0.33130 


CI 


0.43348 


0.43372 


0.43366 


0.43339 


Ar 


0.54298 


0.54444 


0.54433 


0.54395 


Average 


error 


3 X 10-* 


3 X 10"'' 


6 X 10"'' 


Maximum error 


1 X 10"^ 


1 X 10"^ 


-3 X 10"^ 



TABLE IIL Comparison of ionisation potentials obtained 
from AE calculations (E^!^), from calculations using the 
extremely non-local and localised pseudopotentials gener- 
ated within this paper {Eion and Elon respectively), and 
from calculations using the pseudopotentials generated by 
Ovcharenko et a/.— {EZi^) (a.u.). 



